

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
clearvars -except thetaStar K mpow N reb X Y YY Z specif


%% load everything of interest
load('sim_eq_final.mat')
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% model predictions
Yh = Yhat(thetaStar,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
Yh = Yh ./ sum(Yh); % enforce probs add up to 1 (due to Gambit failures)


%% model fit
modelFit=table;
% civil war fit
modelFit.Civil_War=[mean(YY(Y,1));mean(1-Yh(1,:))];
% intervention fit
modelFit.US_Interv=[mean(YY(Y,2));mean(sum(Yh(YY(:,2)==1,:)))];
modelFit.UK_Interv=[mean(YY(Y,3));mean(sum(Yh(YY(:,3)==1,:)))];
modelFit.France_Interv=[mean(YY(Y,4));mean(sum(Yh(YY(:,4)==1,:)))];
modelFit.Russia_Interv=[mean(YY(Y,5));mean(sum(Yh(YY(:,5)==1,:)))];
modelFit.China_Interv=[mean(YY(Y,6));mean(sum(Yh(YY(:,6)==1,:)))];
modelFit.Properties.RowNames={'Data','Model'};
disp(' ')
disp('for Table C1:')
disp(modelFit)


%%
clearvars -except Yh modelFit
save('model_fit.mat')



